{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# 6. Inverse theory methods"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"%plot --format svg -s 800,600"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 6.1 Review of the least-squares fitting\n",
"\n",
"Linear regression that we used for the last project, was essentially a problem of minimizing the following sum of squares of residuals, i.e. the differences between the data $(x_i,y_i)$ and the linear model $y(x)=c_0+c_1 x$ predictions:\n",
"$$\n",
" \\min \\left[ \\sum_{i=1}^m r_i^2 \\right] , \\quad r_i = y_i - (c_0 + c_1 x_i)\n",
"$$\n",
"In matrix form, the problem is to solve the set of linear equations\n",
"$$\n",
" \\vec{y} = \\left( \\begin{array}{c} y_1 \\\\ y_2 \\\\ \\vdots \\\\ y_m\\end{array} \\right) =\n",
" \\left( \\begin{array}{cc} 1 & x_1 \\\\ 1 & x_2 \\\\ \\vdots & \\dots \\\\ 1 & x_m\\end{array} \\right)\n",
" \\left( \\begin{array}{c} c_0 \\\\ c_1 \\end{array} \\right)\n",
"$$\n",
"for $\\{ c_0,c_1\\}$. In octave this can be accomplished by using operator \\ :"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"c =\r\n",
"\r\n",
" -0.55333\r\n",
" 0.98970\r\n",
"\r\n"
]
},
{
"data": {
"image/svg+xml": [
""
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"x_data = [1:10]';\n",
"y_data = [ 0.2 1.0 2.6 3.6 4.9 5.3 6.5 7.8 8.0 9.0 ]';\n",
"A = [ ones(size(x_data)) x_data ];\n",
"c = A \\ y_data\n",
"\n",
"y_lr = A*c;\n",
"plot(x_data,y_data,'r.;data;',\n",
" x_data,y_lr,'g-;linear regression;');\n",
"ylabel('y');\n",
"xlabel('x');\n",
"legend(\"location\",'southeast');\n",
"box on;\n",
"grid on;\n",
"title('Linear regression');"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"There is another way to look at the solution of this overdetermined problem\n",
"$$\n",
" \\vec{y} = {A} \\vec{c}\n",
"$$\n",
"Any full rank matrix $A ∈ R^{m×n}$ , with $m \\geq n$, admits a unique $QR$ factorization\n",
"$A = QR$. It is possible to prove that $A = \\tilde{Q}\\tilde{R}$ where $\\tilde{Q}=Q(1:m.1:n)$ and $\\tilde{R}=R(1:n,1:n)$ are the submatrices as showni below (figure from A.Quarteroni, F.Saleri and P.Gervasio, Scientific Computing with\n",
"MATLAB and Octave, 3rd ed, Springer).\n",
"
\n",
"\n",
"
\n",
"$\\tilde{Q}$ has orthonormal column vectors, while $\\tilde{R}$ is a non-singular upper triangular matrix. Therefore, the unique solution is given by\n",
"$$\n",
" \\vec{c}^* = (\\tilde{R})^{-1} (\\tilde{Q})^T \\vec{y}\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"ans =\n",
"\n",
" 1.00000 0.00000\n",
" 0.00000 1.00000\n",
"\n",
"ans =\n",
"\n",
" 0 0\n",
" 0 0\n",
"\n",
"ans =\n",
"\n",
" -0.55333 -0.55333\n",
" 0.98970 0.98970\n",
"\n"
]
}
],
"source": [
"### another view of the same linear regression problem: QR factorization\n",
"### Again, solving A c = y for c by using A=QR\n",
"[m,n] = size(A);\n",
"### only works for an OVERdetermined system, m>n (here m=10, n=2)\n",
"[Q,R] = qr(A);\n",
"## \\tilde{Q} is an (mxn) subset of Q(mxm), \\tilde{R} is an (nxn) subset of R(mxn)\n",
"Qt = Q(:,1:n); Rt = R(1:n,:);\n",
"## \\tilde{Q} is an orthogonal matrix, so Qt`*Qt should be an identity:\n",
"Qt'*Qt\n",
"## R is already an upper-triungular matrix, so this should be identically zero:\n",
"Rt - triu(Rt)\n",
"## and this yields exactly the same result as before\n",
"[c (Rt \\ (Qt'*y_data)) ]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"But what if the task is not to minimize the (vertical) differences between $y_i$ and $c_0 + c_1 x_i$, but instead the sum of the distances between the data points and the line of fit? Another formulation of the problem is the \"minimum length of water pipes\" of an example distributed with eXtrema, see https://www.physics.brocku.ca/Labs/extrema/.\n",
"\n",
"On a plane, a straight line can be represented by the equation\n",
"$$\n",
" c_0 + c_1 x + c_2 y = 0, \\quad c_1^2 + c_2^2 = 1\n",
"$$\n",
"where $(c_1,c_2)$ is the unit vector perpendicular to the line. For points $(x_i,y_i)$ not on the line, the above equation represents the distances of the points from the line,\n",
"$$\n",
" r_i = c_0 + c_1 x_i + c_2 y_i \n",
"$$\n",
"and the problem becomes the minimization of the total sum of such distances:\n",
"$$\n",
" \\min ||\\vec{r} ||^2 = \\min \\left[ \\sum_{i=1}^m r_i^2 \\right]\n",
"$$\n",
"subject to\n",
"$$\n",
" \\vec{r} = \\left( \\begin{array}{c} r_1 \\\\ r_2 \\\\ \\vdots \\\\ r_m\\end{array} \\right) =\n",
" \\left( \\begin{array}{ccc} 1 & x_1 & y_1 \\\\ 1 & x_2 & y_2 \\\\ \\vdots & \\dots & \\vdots\\\\ 1 & x_m & y_m\\end{array} \\right)\n",
" \\left( \\begin{array}{c} c_0 \\\\ c_1 \\\\ c_2\\end{array} \\right) = A \\vec{c}\n",
" \\quad \\mbox{and} \\quad\n",
" c_1^2+c_2^2=1\n",
"$$\n",
"\n",
"Commands \n",
"qr(), \n",
"triu(), and\n",
"svd()\n",
"will come handy."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"ans =\r\n",
"\r\n",
" -0.41624 0.70565 -0.70856 1.00000 1.00000 1.00000\r\n",
"\r\n"
]
},
{
"data": {
"image/svg+xml": [
""
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# this function solves the constrained least squares problem r = A c\n",
"A = [ ones(size(x_data)) x_data y_data ];\n",
"## m data, n parameters of the fit\n",
"[m,n] = size(A);\n",
"## on a two-dimensional plane, line normal only has two \n",
"## components, (c_1,c_2), so only two parameters of the fit\n",
"dim = 2; \n",
"if n < dim+1, error ('not enough parameters'); end;\n",
"if m < dim, error ('not enough data points'); end;\n",
"\n",
"## QR factorization\n",
"[Q,R] = qr(A);\n",
"## Singular value decomposition: V^T S* U = pseudoinverse of A\n",
"[U,S,V] = svd(R(n-dim+1:m,n-dim+1:n)); \n",
"cc = V(:,dim); ## only dim constraints are involved\n",
"cc_0 = -R(1:n-dim,1:n-dim) \\ (R(1:n-dim,n-dim+1:n)*cc);\n",
"\n",
"[ [cc_0 cc'] (sqrt(cc(1)^2+cc(2)^2)) cc'*cc norm(cc)]\n",
"\n",
"y_ls = -(cc_0+cc(1)*x_data)/cc(2);\n",
"plot(x_data,y_data,'r.;data;',\n",
" x_data,y_lr,'g-;linear regression;',\n",
" x_data,y_ls,'b-;constrained least squares;',\n",
" x_data,100*(y_lr - y_ls),'c-;100 x difference;');\n",
"ylabel('y');\n",
"xlabel('x');\n",
"legend(\"location\",'east');\n",
"box on;\n",
"grid on;\n",
"title(sprintf('Least-squares fitting: c = (%.3f, %.3f, %.3f)',cc_0,cc(1),cc(2)));"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 6.2 Introduction to inverse problems"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"scrolled": true
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"08_Sternin-reprint.pdf\t60C.dat expo_maketest.m\t pmma.eps test.dat\r\n",
"30C.dat\t\t\t70C.dat expo_maketest.sci\t pmma.png\r\n",
"40C.dat\t\t\t80C.dat expo.sci\t\t pmma.ps\r\n",
"50C.dat\t\t\texpo.m\t inverse-theory-methods.pdf SS97.pdf\r\n"
]
}
],
"source": [
"%bash\n",
"ls /work/5P10/Inverse/"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Inverse Theory presentation (pdf)"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
""
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"image/svg+xml": [
""
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"%% expo_maketest.m\n",
"% make f(t) + noise, for testing TR analysis of \n",
"% a distribution of relaxation rates, g(r) \n",
"%\n",
"% Completed: Dec.2005, Edward Sternin \n",
"% Contributions: Hartmut Schaefer\n",
"% Revisions: 2018.02 ES converted to matlab/octave from SciLab\n",
"%\n",
"clear;\n",
"\n",
"noise=0.05; % noise as fraction of the max(f(t))\n",
"datafile=sprintf(\"noisy-%.2f.dat\",noise);\n",
"\n",
"t_steps=1000; % t-domain signal has this many points,\n",
"t_min=1; % from this minimum,\n",
"t_max=1e5; % to this maximum\n",
"r_steps=500; % r-grid has this many points,\n",
"r_min=1/t_max; % from this minimum,\n",
"r_max=1/t_min; % to this maximum\n",
"\n",
"% Use an array of Gaussian peaks that make up the true g(r):\n",
"%peaks = [[0.5 0.20 0.04];[0.65 0.50 0.03];[0.75 0.80 0.03]]; % [amplitude center width]\n",
"peaks = [[0.5 0.20 0.04];[0.65 0.50 0.03]]; % [amplitude center width]\n",
"%peaks = [[0.6 0.35 0.1];[0.1 0.45 0.15]]; % [amplitude center width]\n",
"\n",
"% create a vector of r values, logarithmically spaced\n",
"r_scale = log(r_max)-log(r_min);\n",
"r = exp([log(r_min):r_scale/(r_steps-1):log(r_max)])'; % use a column vector\n",
"\n",
"% build up g(r) by adding Gaussian contributions from all peaks\n",
"g = zeros(length(r),1);\n",
"for k=1:length(peaks(:,1))\n",
" r0 = log(r_min) + r_scale*peaks(k,2);\n",
" dr = r_scale*peaks(k,3);\n",
" g += peaks(k,1) * exp (- (log(r) - r0).^2 / (2 * dr^2));\n",
"end\n",
"\n",
"f1=figure(1); clf(1);\n",
"semilogx(r,g,'r-','LineWidth',2);\n",
"ylabel(\"g(r)\");\n",
"xlabel(\"r\");\n",
"title(\"Test distribution of relaxation rates\");\n",
"FS = findall(f1,'-property','FontSize');\n",
"set(FS,'FontSize',14);\n",
"%print -dpng g.png\n",
"\n",
"% create a vector of t values, linearly spaced; use column vectors for f(t)\n",
"t = [t_min:(t_max-t_min)/(t_steps-1):t_max]';\n",
"\n",
"% generate time-domain data\n",
"f = (1/(length(t)-1)*g'*exp(-r*t'))';\n",
"\n",
"% normalize and add random noise\n",
"f /= max(f);\n",
"f += noise*(-0.5+rand(length(f),1));\n",
"\n",
"f2=figure(2); clf(2);\n",
"plot(t,f,'g-');\n",
"%semilogy(t,f,'g-');\n",
"ylabel(\"f(t)\");\n",
"xlabel(\"t\");\n",
"title(sprintf(\"f(t) + %.1f%% random noise, saved in %s\",100*noise,datafile));\n",
"FS = findall(f2,'-property','FontSize');\n",
"set(FS,'FontSize',14);\n",
"%print -dpng f.png\n",
"\n",
"%system(sprintf(\"rm -rf %s\",datafile));\n",
"f_of_t=[t f];\n",
"save(\"-ascii\",datafile,\"f_of_t\");\n",
"\n",
"g_of_r=[r g];\n",
"save(\"-ascii\",\"true_g.dat\",\"g_of_r\");\n"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
""
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"image/svg+xml": [
"
![\"[Fig5.9]\"](\"QRfactorization.png\")
\n",
"